Read in Data

Glasser regions and assignments

#Glasser region and label names for the frontal lobe
glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")

glasser.plotting <- glasser.regions
glasser.plotting$cortex <- "cortex"
glasser.plotting$cortex[glasser.plotting$orig_parcelname %in% glasser.snr.exclude$orig_parcelname] <- "exclude"
glasser.plotting <- glasser.plotting %>% select(label, cortex)
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "rh_???", cortex = "exclude"))
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "lh_???", cortex = "exclude"))
glasser.plotting <- glasser.plotting %>% filter(label != "lh_L_10pp")

Depths

depth.list <- c("depth_1", "depth_2", "depth_3", "depth_4", "depth_5", "depth_6", "depth_7")
ordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <- c("#004A38FF", "#006458FF", "#0093B0FF", "#8AABD5FF", "#C0BFE3FF", "#E0D5ECFF", "#F2E7F4FF")

S-A Axis

#S-A axis or saxis as nicknamed by Dan Margulies on a trip to gingerworld
S.A.axis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = "orig_parcelname")

Neurosynth Term Scores

neurosynth.terms <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Neurosynth_terms/atl-glasser360_neurosynth.csv") %>% dplyr::select(-regionName) %>% dplyr::select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of term names
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rh
neurosynth.terms[,1:123] <- scale(neurosynth.terms[,1:123])

Final participant list

participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv")
participants <- participants %>% mutate(subses = sprintf("%s_%s", subject_id, session_id))

Depth-dependent R1 in HCP-MMP regions

#R1 by region matrices at 7 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.R

GAM outputs: developmental effects

setwd("/Volumes/Hera/Projects/corticalmyelin_development/gam_outputs/developmental_effects/") #output from /gam_models/fitgams_glasserregions_bydepth.R

files <- list.files(getwd(), pattern = "age.RDS") 

#read in files and assign to variables
for(i in 1:length(files)){
  
  Rfilename <- gsub(".RDS", "", files[i])
  
  x <- readRDS(files[i]) 
  assign(Rfilename, x) 
  }

#Combine all outputs stats. A list of lists! #inception #spinning top
#7 depths included, each depth contains 4 dfs (stats, fitted, smooths, derivatives)
gam.results.alldepths <- list(depth1_gamstatistics_age, depth2_gamstatistics_age, depth3_gamstatistics_age, depth4_gamstatistics_age, depth5_gamstatistics_age, depth6_gamstatistics_age, depth7_gamstatistics_age)
names(gam.results.alldepths) <- list("depth_1", "depth_2", "depth_3", "depth_4", "depth_5", "depth_6", "depth_7")

#Extract the 4 types of results and format
gam.statistics.alldepths <- lapply(gam.results.alldepths, '[[', "gam.statistics.df" )
gam.statistics.allregions <- gam.statistics.alldepths #whole brain data
gam.fittedvalues.alldepths <-  lapply(gam.results.alldepths, '[[', "gam.fittedvalues.df" )
gam.smoothestimates.alldepths <- lapply(gam.results.alldepths, '[[', "gam.smoothestimates.df" )
gam.smoothestimates.allregions<- lapply(gam.results.alldepths, '[[', "gam.smoothestimates.df" ) #wholebrain data
gam.derivatives.alldepths <- lapply(gam.results.alldepths, '[[', "gam.derivatives.df" )

format_depth_dfs <- function(depth){
  depth <- depth[depth$orig_parcelname %in% glasser.frontal$orig_parcelname,] #get results for frontal lobe only
  depth <-  depth[!(depth$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
}

gam.statistics.alldepths <- lapply(gam.statistics.alldepths, function(depth){
  format_depth_dfs(depth)
})
gam.fittedvalues.alldepths <- lapply(gam.fittedvalues.alldepths, function(depth){
  format_depth_dfs(depth)
})
gam.smoothestimates.alldepths <- lapply(gam.smoothestimates.alldepths, function(depth){
  format_depth_dfs(depth)
})
gam.derivatives.alldepths <- lapply(gam.derivatives.alldepths, function(depth){
  format_depth_dfs(depth)
})

Myelin Matures Heterochronously and Hierarchically Across Prefrontal Cortical Depths

Frontal lobe developmental trajectory at each depth

#Long formatted R1 data for all scans + frontal lobe regions at each depth
myelin.glasser.7T.long <- lapply(myelin.glasser.7T, function(depth){
  cols_to_pivot <- names(depth)[grepl("ROI", names(depth))]
  depth.long <- depth %>% pivot_longer(cols = cols_to_pivot, names_to = "orig_parcelname", values_to = "R1")
  depth.long <- merge(depth.long, glasser.frontal, by = "orig_parcelname")
  depth.long <- depth.long[!(depth.long$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
  depth.long <- depth.long %>% mutate(subses = sprintf("%s_%s", subject_id, session_id))
  return(depth.long)
})

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(cols_to_pivot)
## 
##   # Now:
##   data %>% select(all_of(cols_to_pivot))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Calculate average frontal lobe R1 for each sub/ses
myelin.frontallobe.7T <- lapply(myelin.glasser.7T.long, function(depth){
    depth <- depth %>% group_by(subses) %>% do(mean.R1 = mean(.$R1)) %>% unnest(cols = c( "mean.R1"))
    depth <- merge(depth, participants, by = "subses")
    depth$subject_id <- as.factor(depth$subject_id)
    depth$sex <- as.factor(depth$sex)
    return(depth)
})

#Function to fit a frontal lobe gam for a given depth and return model summary
gam.summary <- function(depth){
  
  depth.data <- depth
  
  model <- gam(mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re"), method = c("REML"), data = depth.data)
  
  print(summary(model))
}

#Function to fit a frontal lobe gam for a given depth and plot the developmental trajectory
plot.trajectory <- function(depth, display_color, y_lims, y_ticks){
  
  depth.data <- depth
  
  #First fit the gam and get fitted values and derivatives 
  depth.model <- gam.statistics.smooths(input.df = depth.data, region = "mean.R1", smooth_var = "age", id_var = "subject_id", covariates = "NA", random_intercepts = T, random_slopes = F, knots = 4, set_fx = F)

  #Age spline plot with participant-level data
  trajectory.plot <- ggplot() +
    geom_point(data = depth.data, aes(x = age, y = mean.R1, group = subject_id), color = "black", size = .4) +
    geom_line(data = depth.data, aes(x = age, y = mean.R1, group = subject_id), linewidth = 0.3, color = "gray75") +
    geom_line(data = depth.model$gam.fittedvalues, aes(x = age, y = fitted), color = display_color, linewidth = 1.5) +
    labs(x="\nAge", y=sprintf("R1\n")) +
    theme_classic() +
    scale_y_continuous(limits = y_lims, breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(trajectory.plot)
}

#Function to fit a frontal lobe gam for a given depth and return the fitted values df
depth.fittedvalues <- function(depth){
  
  depth.data <- depth
  
  #First fit the gam and get fitted values
  depth.model <- gam.statistics.smooths(input.df = depth.data, region = "mean.R1", smooth_var = "age", id_var = "subject_id", covariates = "NA", random_intercepts = T, random_slopes = F, knots = 4, set_fx = F)
  
  depth.fitted <- depth.model$gam.fittedvalues
  return(depth.fitted)
}

Depth-speficic plots

Depth 1

gam.summary(myelin.frontallobe.7T$depth_1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5145777  0.0008383   613.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf  Ref.df      F  p-value    
## s(age)         1.236   1.357 24.720 4.30e-07 ***
## s(subject_id) 67.102 139.000  1.049 3.25e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.476   Deviance explained = 64.4%
## -REML = -663.39  Scale est. = 7.0085e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_1, display_color = "#EBE0F0FF", y_lims = c(0.48, 0.555), y_ticks  = c(0.49, 0.52, 0.55))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth1_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 2

gam.summary(myelin.frontallobe.7T$depth_2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5261190  0.0008012   656.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf  Ref.df      F  p-value    
## s(age)         1.943   2.216 31.910  < 2e-16 ***
## s(subject_id) 62.903 139.000  0.916 1.69e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.525   Deviance explained = 66.9%
## -REML = -670.04  Scale est. = 6.7933e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_2, display_color = depth_colorbar[6], y_lims = c(0.495, 0.56), y_ticks = c(0.50, 0.53, 0.56))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth2_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 3

gam.summary(myelin.frontallobe.7T$depth_3)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5364280  0.0007889     680   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf  Ref.df      F  p-value    
## s(age)         2.199   2.473 42.504  < 2e-16 ***
## s(subject_id) 60.074 139.000  0.833 5.04e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.57   Deviance explained = 69.5%
## -REML = -671.36  Scale est. = 6.8544e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_3, display_color = depth_colorbar[5], y_lims = c(0.505, 0.57), y_ticks = c(0.51, 0.54, 0.57))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth3_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 4

gam.summary(myelin.frontallobe.7T$depth_4)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.546799   0.000796   686.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf  Ref.df      F  p-value    
## s(age)         2.325   2.591 51.691  < 2e-16 ***
## s(subject_id) 57.575 139.000  0.768 0.000121 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =    0.6   Deviance explained = 71.2%
## -REML = -667.8  Scale est. = 7.226e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_4, display_color = depth_colorbar[4], y_lims = c(0.51, 0.585), y_ticks = c(0.52, 0.55, 0.58))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth4_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 5

gam.summary(myelin.frontallobe.7T$depth_5)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5596640  0.0008176   684.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf  Ref.df      F  p-value    
## s(age)         2.39   2.649 58.508  < 2e-16 ***
## s(subject_id) 55.42 139.000  0.716 0.000238 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.618   Deviance explained = 72.1%
## -REML = -660.76  Scale est. = 7.8505e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_5, display_color = depth_colorbar[3], y_lims = c(0.525, 0.6), y_ticks = c(0.53, 0.56, 0.59))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth5_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 6

gam.summary(myelin.frontallobe.7T$depth_6)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5779681  0.0008437   685.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf  Ref.df      F  p-value    
## s(age)         2.414   2.669 62.514  < 2e-16 ***
## s(subject_id) 54.658 139.000  0.698 0.000304 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.629   Deviance explained = 72.8%
## -REML = -653.6  Scale est. = 8.4454e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_6, display_color = depth_colorbar[2], y_lims = c(0.54, 0.63), y_ticks = c(0.56, 0.59, 0.62))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth6_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

Depth 7

gam.summary(myelin.frontallobe.7T$depth_7)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.6041477  0.0008685   695.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf  Ref.df      F p-value    
## s(age)         2.40   2.654 62.389 < 2e-16 ***
## s(subject_id) 57.61 139.000  0.764 0.00014 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.641   Deviance explained = 74.2%
## -REML = -649.14  Scale est. = 8.5938e-05  n = 215

plot.trajectory(depth = myelin.frontallobe.7T$depth_7, display_color = depth_colorbar[1], y_lims = c(0.565, 0.65), y_ticks = c(0.58, 0.61, 0.64))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth7_trajectoryplot.pdf", device = "pdf", dpi = 500, width = 2.48, height = 1.48)

ps <- c(4.30e-07, 2e-16, 2e-16, 2e-16, 2e-16, 2e-16, 2e-16)
p.adjust(ps, method = c("fdr"))

## [1] 4.300000e-07 2.333333e-16 2.333333e-16 2.333333e-16 2.333333e-16
## [6] 2.333333e-16 2.333333e-16

Multi-depth plots

frontal.fittedvalues.alldepths <- lapply(myelin.frontallobe.7T, function(depth){
  depth.fittedvalues(depth)
})

ggplot() +
#depth 1
geom_point(data = myelin.frontallobe.7T$depth_1, aes(x = age, y = mean.R1, group = subject_id), color = "black", size = .5) +
geom_line(data = myelin.frontallobe.7T$depth_1, aes(x = age, y = mean.R1, group = subject_id), linewidth = 0.3, color = "gray75") +
geom_line(data = frontal.fittedvalues.alldepths$depth_1, aes(x = age, y = fitted), color = depth_colorbar[7], linewidth = 1.8) +
#depth 7
geom_point(data = myelin.frontallobe.7T$depth_7, aes(x = age, y = mean.R1, group = subject_id), color = "black", size = .5) +
geom_line(data = myelin.frontallobe.7T$depth_7, aes(x = age, y = mean.R1, group = subject_id), linewidth = 0.3, color = "gray75") +
geom_line(data = frontal.fittedvalues.alldepths$depth_7, aes(x = age, y = fitted), color = depth_colorbar[1], linewidth = 1.8) +
labs(x="\nAge", y=sprintf("R1\n")) +
theme_classic() +
scale_y_continuous(limits = c(0.48, 0.65)) +
theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2_supplementary/Depth_trajectories_withdata.pdf", device = "pdf", dpi = 500, width = 2.7, height = 4.65)

ggplot() +
#depth 1
geom_line(data = frontal.fittedvalues.alldepths$depth_1, aes(x = age, y = fitted), color = depth_colorbar[7], linewidth = 1.8) +
#depth 2
geom_line(data = frontal.fittedvalues.alldepths$depth_2, aes(x = age, y = fitted), color = depth_colorbar[6], linewidth = 1.8) +
#depth 3
geom_line(data = frontal.fittedvalues.alldepths$depth_3, aes(x = age, y = fitted), color = depth_colorbar[5], linewidth = 1.8) +
#depth 4
geom_line(data = frontal.fittedvalues.alldepths$depth_4, aes(x = age, y = fitted), color = depth_colorbar[4], linewidth = 1.8) +
#depth 5
geom_line(data = frontal.fittedvalues.alldepths$depth_5, aes(x = age, y = fitted), color = depth_colorbar[3], linewidth = 1.8) +
#depth 6
geom_line(data = frontal.fittedvalues.alldepths$depth_6, aes(x = age, y = fitted), color = depth_colorbar[2], linewidth = 1.8) +
#depth 7
geom_line(data = frontal.fittedvalues.alldepths$depth_7, aes(x = age, y = fitted), color = depth_colorbar[1], linewidth = 1.8) +
labs(x="\nAge", y=sprintf("R1\n")) +
theme_classic() +
scale_y_continuous(limits = c(0.48, 0.65)) +
theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2_supplementary/Depth_trajectories_allfits.pdf", device = "pdf", dpi = 500, width = 2.7, height = 4.65)

Number of significant regional effects at each depth

#Function to calculate the number of significant age effects post-FDR
ageeffect.results <- function(depth){
  ageeffects <- depth
  n_roi <- nrow(ageeffects)
 
  #significant smooth terms
  aageeffects <- ageeffects %>% mutate(significant = p.adjust(ageeffects$GAM.smooth.pvalue, method = c("fdr")) < 0.05)
  sig_age <- aageeffects %>% filter(significant == TRUE) %>% nrow()
  sig_age_percent <- round((sig_age/n_roi)*100, 2)
                         
  sprintf("%s percent of regions show a positive significant developmental change in R1", sig_age_percent)
}

lapply(gam.statistics.alldepths, function(depth){
  ageeffect.results(depth)
})

## $depth_1
## [1] "85.93 percent of regions show a positive significant developmental change in R1"
## 
## $depth_2
## [1] "92.59 percent of regions show a positive significant developmental change in R1"
## 
## $depth_3
## [1] "92.59 percent of regions show a positive significant developmental change in R1"
## 
## $depth_4
## [1] "95.56 percent of regions show a positive significant developmental change in R1"
## 
## $depth_5
## [1] "97.78 percent of regions show a positive significant developmental change in R1"
## 
## $depth_6
## [1] "98.52 percent of regions show a positive significant developmental change in R1"
## 
## $depth_7
## [1] "98.52 percent of regions show a positive significant developmental change in R1"

Regional derivatives at each depth

gam.derivatives.alldepths.long <- do.call(rbind, gam.derivatives.alldepths)
gam.derivatives.alldepths.long <- gam.derivatives.alldepths.long %>% mutate(depth = as.factor(substr(row.names(gam.derivatives.alldepths.long), 1, 7))) #add depth as a column
gam.derivatives.alldepths.long <- merge(gam.derivatives.alldepths.long, glasser.regions, sort = F)

#Calculate the average first derivative in each region at each depth
regional.rate.alldepths <- gam.derivatives.alldepths.long %>% group_by(label, depth) %>% do(rate = mean(.$derivative)) %>% unnest(cols = "rate")

#Calculate the average across-region first derivative at each depth
frontallobe.rate.alldepths <- regional.rate.alldepths %>% group_by(depth) %>% do(mean.rate = mean(.$rate)) %>% unnest(cols = "mean.rate")
frontallobe.rate.alldepths$depth <- factor(frontallobe.rate.alldepths$depth, ordered = T, levels = ordered_depths)

Depth-wise derivative plot

frontallobe.rate.alldepths %>% 
  ggplot(aes(x = mean.rate, y = depth, fill = depth)) +
  geom_point(shape = 23, size = 5.5, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  xlab("\nMean Derivative") +
  ylab("Cortical Depth\n") +
  scale_x_continuous(breaks = c(0.0010, 0.0014, 0.0018), limits = c(0.00085, 0.0019)) +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depthplot_derivative.pdf", device = "pdf", dpi = 500, width = 2.1, height = 4.0)

Depth-wise derivative brain plots

for(this.depth in c(unique(regional.rate.alldepths$depth))){
  
  regional.rate.depth <- regional.rate.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  derivative.plot <- ggplot() + 
  geom_brain(data = regional.rate.depth, atlas = glasser, mapping = aes(fill = rate), colour=I("gray50"), size=I(.08)) + theme_classic() + theme(legend.position = "none") + 
  paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.0007, 0.0022), oob = squish, na.value = "white") +  
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 
  
  print(derivative.plot)
  
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/%s_rate_corticalmap.pdf", this.depth), device = "pdf", dpi = 500, width = 3.95, height = 2)
}

Depth correlation matrix

regional.rate.alldepths.wide <- regional.rate.alldepths %>% pivot_wider(id_cols = label, names_from = depth, values_from = rate) %>% select(-label)
regional.rate.alldepths.wide <- regional.rate.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1) #order for plot
depth.derivative.cors <- cor(regional.rate.alldepths.wide) # compute correlation matrix

ggcorrplot(corr = depth.derivative.cors, method = "square", type = "upper", show.diag = T, lab = F) +
scale_fill_gradientn(colors = c("white", "#CCDCED", "#BFD2E2", "#7B9FC6"), limits = c(0.7, 1))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth_correlationmatrix_derivative.pdf", device = "pdf", dpi = 500, width = 2.85, height = 2)

Coefficient of variation

rate.cov.bydepth <- data.frame(depth = factor(), cov = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.rate.alldepths %>% filter(depth == this.depth)
  cov <- cv(depth.data$rate)
  cov <- data.frame(depth = this.depth, cov = cov)
  rate.cov.bydepth <- rbind(rate.cov.bydepth, cov)
}
rate.cov.bydepth

##     depth       cov
## 1 depth_1 0.5074685
## 2 depth_2 0.4501392
## 3 depth_3 0.4168052
## 4 depth_4 0.3802468
## 5 depth_5 0.3417924
## 6 depth_6 0.3086435
## 7 depth_7 0.2845120

S-A axis correlation

rate.SAaxis.bydepth <- data.frame(depth = factor(), SA.corr = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.rate.alldepths %>% filter(depth == this.depth)
  depth.data <- merge(depth.data, S.A.axis)
  corr <- cor(depth.data$SA.axis, depth.data$rate, method = c("spearman"), use = "complete.obs")
  corr <- data.frame(depth = this.depth, SA.corr = corr)
  rate.SAaxis.bydepth <- rbind(rate.SAaxis.bydepth, corr)
}
rate.SAaxis.bydepth

##     depth    SA.corr
## 1 depth_1 -0.5096381
## 2 depth_2 -0.4888304
## 3 depth_3 -0.4633743
## 4 depth_4 -0.4102722
## 5 depth_5 -0.3581602
## 6 depth_6 -0.2967906
## 7 depth_7 -0.2408838

rate.SAaxis.bydepth$depth <- factor(rate.SAaxis.bydepth$depth, ordered = T, levels = ordered_depths)
ggplot(rate.SAaxis.bydepth, aes(x = SA.corr, y = depth,  fill = depth)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  xlab("\nS-A Axis Correlation") +
  ylab("Cortical Depth\n") +
  scale_x_continuous(breaks = c(-0.5, -0.25, 0), limits = c(-0.52, 0)) +
  scale_fill_manual(values = depth_colorbar) +
  theme(axis.text.y = element_text(angle = 45, hjust = 1)) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Derivative_SAaxis_plot.pdf", device = "pdf", dpi = 500, width = 2.05, height = 1.55)

Regional age of maturation at each depth

gam.statistics.alldepths.long <- do.call(rbind, gam.statistics.alldepths)
gam.statistics.alldepths.long <- gam.statistics.alldepths.long %>% mutate(depth = as.factor(substr(row.names(gam.statistics.alldepths.long), 1, 7))) #add depth as a column
regional.maturation.alldepths <- gam.statistics.alldepths.long
regional.maturation.alldepths <- merge(regional.maturation.alldepths, glasser.regions, sort = F)

#Calculate the average across-region age of maturation at each depth
frontallobe.maturation.alldepths <- regional.maturation.alldepths %>% group_by(depth) %>% do(mean.age = mean(.$smooth.increase.offset, na.rm = T)) %>% unnest(cols = "mean.age")
frontallobe.maturation.alldepths$depth <- factor(frontallobe.maturation.alldepths$depth, ordered = T, levels = ordered_depths)

Depth-wise age of maturation plot

frontallobe.maturation.alldepths %>% 
  ggplot(aes(x = mean.age, y = depth, fill = depth)) +
  geom_point(shape = 23, size = 5.5, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  xlab("\nAge of Maturation") +
  ylab("Cortical Depth\n") +
  scale_x_continuous(limits = c(25.95, 30.5)) +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depthplot_maturation.pdf", device = "pdf", dpi = 500, width = 2.1, height = 4.0)

Depth-wise age of maturation brain plots

for(this.depth in c(unique(regional.maturation.alldepths$depth))){
  
  regional.mat.depth <- regional.maturation.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  derivative.plot <- ggplot() + 
  geom_brain(data = regional.mat.depth, atlas = glasser, mapping = aes(fill = smooth.increase.offset), colour=I("gray50"), size=I(.08)) + theme_classic() + theme(legend.position = "none") + 
  paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = 1, limits = c(25.5, 33.5), oob = squish, na.value = "white") +  
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 
  
  print(derivative.plot)
  
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/%s_maturation_corticalmap.pdf", this.depth), device = "pdf", dpi = 500, width = 3.95, height = 2)
}

Depth correlation matrix

regional.maturation.alldepths.wide <- regional.maturation.alldepths %>% pivot_wider(id_cols = label, names_from = depth, values_from = smooth.last.change) %>% select(-label)
regional.maturation.alldepths.wide <- regional.maturation.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1) #order for plot
depth.maturation.cors <- cor(regional.maturation.alldepths.wide, use = "complete.obs") # compute correlation matrix

ggcorrplot(corr = depth.maturation.cors, method = "square", type = "upper", show.diag = T, lab = F) +
scale_fill_gradientn(colors = c("white", "#CCDCED", "#BFD2E2", "#7B9FC6"), limits = c(0.3, 1))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Depth_correlationmatrix_maturation.pdf", device = "pdf", dpi = 500, width = 2.85, height = 2)

Coefficient of variation

maturation.cov.bydepth <- data.frame(depth = factor(), cov = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.maturation.alldepths %>% filter(depth == this.depth)
  cov <- cv(depth.data$smooth.increase.offset, na.rm = T)
  cov <- data.frame(depth = this.depth, cov = cov)
  maturation.cov.bydepth <- rbind(maturation.cov.bydepth, cov)
}
maturation.cov.bydepth

##     depth       cov
## 1 depth_1 0.1297405
## 2 depth_2 0.1520745
## 3 depth_3 0.1417636
## 4 depth_4 0.1372049
## 5 depth_5 0.1365760
## 6 depth_6 0.1350223
## 7 depth_7 0.1286170

S-A axis correlation

maturation.SAaxis.bydepth <- data.frame(depth = factor(), SA.corr = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.maturation.alldepths %>% filter(depth == this.depth)
  depth.data <- merge(depth.data, S.A.axis)
  corr <- cor(depth.data$SA.axis, depth.data$smooth.increase.offset, method = c("spearman"), use = "complete.obs")
  corr <- data.frame(depth = this.depth, SA.corr = corr)
  maturation.SAaxis.bydepth <- rbind(maturation.SAaxis.bydepth, corr)
}
maturation.SAaxis.bydepth

##     depth    SA.corr
## 1 depth_1 0.24376118
## 2 depth_2 0.21074195
## 3 depth_3 0.08614598
## 4 depth_4 0.02608714
## 5 depth_5 0.02949916
## 6 depth_6 0.04376520
## 7 depth_7 0.04846727

maturation.SAaxis.bydepth$depth <- factor(maturation.SAaxis.bydepth$depth, ordered = T, levels = ordered_depths)
ggplot(maturation.SAaxis.bydepth, aes(x = SA.corr, y = depth,  fill = depth)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  xlab("\nS-A Axis Correlation") +
  ylab("Cortical Depth\n") +
  scale_x_continuous(breaks = c(0, .1, .2, .3), limits = c(0, 0.3)) +
  scale_fill_manual(values = depth_colorbar) +
  theme(axis.text.y = element_text(angle = 45, hjust = 1)) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure2/Maturation_SAaxis_plot.pdf", device = "pdf", dpi = 500, width = 2.05, height = 1.55)

Depth-dependent Differences in Myelin Maturational Profiles Vary Across Prefrontal Cortical Zones

Depth-dependent trajectories in exemplar frontal regions

gam.smoothestimates.alldepths.long <- do.call(rbind, gam.smoothestimates.alldepths)
gam.smoothestimates.alldepths.long <- gam.smoothestimates.alldepths.long %>% mutate(depth = as.factor(substr(row.names(gam.smoothestimates.alldepths.long), 1, 7))) #add depth as a column
gam.smoothestimates.alldepths.long$depth <- factor(gam.smoothestimates.alldepths.long$depth, ordered = T, levels = ordered_depths)

#Function to plot age smooth functions in exemplar regions
plot.depth.smooths <- function(region, y_ticks){

  smooth.plot <- gam.smoothestimates.alldepths.long %>% filter(orig_parcelname == region) %>%
  ggplot(., aes(x = age, y = est, group = depth, color = depth)) +
  geom_line(linewidth = .4) +
  theme_classic() +
  xlab("\nAge (years)") +
  ylab("R1 Trajectory (zero-centered)\n") +
  scale_y_continuous(breaks = y_ticks) +
  scale_color_manual(values = depth_colorbar) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(smooth.plot)
}

Primary Motor (4)

glasser.regions %>% filter(orig_parcelname == "L_4_ROI") %>%
ggseg(.data = ., atlas = "glasser", mapping=aes(fill = orig_parcelname, colour=I("gray25"), size=I(.07)), position = c("stacked"), hemisphere = "left", view = "lateral") + 
  theme_void() + 
  labs(fill="") +
  scale_fill_manual(values = c("#395F94"), na.value = "white") +
  theme(legend.position = "none")

## merging atlas and data by 'label'

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Area4_corticalmap.pdf", device = "pdf", dpi = 500, width = .75, height = .55)

plot.depth.smooths(region = "L_4_ROI", y_ticks = c(-0.04, -0.02, 0, 0.02))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Area4_depthsmooths.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.8)

Dorsolateral PFC (46)

glasser.regions %>% filter(orig_parcelname == "L_46_ROI") %>%
ggseg(.data = ., atlas = "glasser", mapping=aes(fill = orig_parcelname, colour=I("gray25"), size=I(.07)), position = c("stacked"), hemisphere = "left", view = "lateral") + 
  theme_void() + 
  labs(fill="") +
  scale_fill_manual(values = c("#395F94"), na.value = "white") +
  theme(legend.position = "none")

## merging atlas and data by 'label'

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Area46_corticalmap.pdf", device = "pdf", dpi = 500, width = .75, height = .55)

plot.depth.smooths(region = "L_46_ROI", y_ticks = c(-0.02, -0.01, 0, 0.01))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Area46_depthsmooths.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.8)

Anterior Cingulate (a24)

glasser.regions %>% filter(orig_parcelname == "L_a24_ROI") %>%
ggseg(.data = ., atlas = "glasser", mapping=aes(fill = orig_parcelname, colour=I("gray25"), size=I(.07)), position = c("stacked"), hemisphere = "left", view = "medial") + 
  theme_void() + 
  labs(fill="") +
  scale_fill_manual(values = c("#395F94"), na.value = "white") +
  theme(legend.position = "none")

## merging atlas and data by 'label'

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Areaa24_corticalmap.pdf", device = "pdf", dpi = 500, width = .75, height = .55)

plot.depth.smooths(region = "L_a24_ROI", y_ticks = c(-0.01, 0, 0.01))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Areaa24_depthsmooths.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.8)

Developmental trajectory curvatures

#A function to compute the curvature of the estimated trajectory
curvature <- function(df, x,  y){
age.spline <- with(df, smooth.spline(x, y, df = 10))
first.deriv <- with(df, predict(age.spline, x = x, deriv = 1))
second.deriv <- with(df, predict(age.spline, x = x, deriv = 2))

k <- (second.deriv$y / ((1 + (first.deriv$y^2))^(3/2)))
k <- abs(k)
return(k)}

#Calculate average curvature of R1 growth trajectories for every region at every cortical depth
trajectory.curvatures <- map_dfr(unique(gam.statistics.alldepths.long$orig_parcelname), function(r){
  region.smooths <- gam.smoothestimates.alldepths.long %>% filter(orig_parcelname == r) #smooth fits for this region at all depths
  depth.curvature <- region.smooths %>% group_by(depth) %>% do(curvature = mean(curvature(df = ., x = .$age , y = .$est))) %>% unnest(cols = c("curvature")) #calculate curvature of the smooth function
  depth.curvature$orig_parcelname <- r
  return(depth.curvature)
})

trajectory.curvatures <- merge(trajectory.curvatures, glasser.regions)
trajectory.curvatures.wide <- trajectory.curvatures %>% pivot_wider(id_cols = c("label"), names_from = "depth", values_from = "curvature") 
labels <- trajectory.curvatures.wide$label
trajectory.curvatures.wide <- trajectory.curvatures.wide %>% dplyr::select(-label)

Assess cluster-ability

library(factoextra)
library("NbClust")
#Assess clustering tendency before applying clustering!
##To assess the clustering tendency, the Hopkins’ statistic is used. It measures the probability that data are generated by a uniform data distribution 
##If the value of Hopkins statistic is close to 1 (far above 0.5), then we can conclude that the dataset is significantly clusterable.

(trajectory.curvatures.wide %>% get_clust_tendency(n = nrow(trajectory.curvatures.wide)-1, seed = 123))$hopkins_stat #hopkins stat + dissimilarity matrix image confirm data is clusterable

## [1] 0.8888908

Determine the optimal number of clusters from a variety of approaches

set.seed(123)

#Use NbClust to determine the best number of clusters to use. This package compares 30 different indexes for determining cluster number and summarizes the best option across all 30

pdf(file = NULL) #to suppress unnescessary plotting

res.nbclust <- trajectory.curvatures.wide %>% scale() %>%
  NbClust(distance = "euclidean",
          min.nc = 2, max.nc = 8, 
          method = "kmeans", index ="all")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
##

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 5 proposed 2 as the best number of clusters 
## * 10 proposed 3 as the best number of clusters 
## * 1 proposed 4 as the best number of clusters 
## * 2 proposed 6 as the best number of clusters 
## * 1 proposed 7 as the best number of clusters 
## * 4 proposed 8 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************

dev.off()

## quartz_off_screen 
##                 2

#And confirm with a scree plot 
set.seed(123)
n_clusters <- 8
wss <- numeric(n_clusters)

#loop over 1 to n possible clusters
for (i in 1:n_clusters) {
  km.out <- kmeans(trajectory.curvatures.wide, centers = i)
  wss[i] <- km.out$tot.withinss
}
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)

ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_line(linewidth = .8, color = depth_colorbar[4]) +
geom_point(size = 1) +
  theme_classic() +
  xlab("\nCluster Number") +
  ylab("Within-cluster SS\n") +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6,7, 8)) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Scree_plot.pdf", device = "pdf", dpi = 500, width = 1.32, height = 1.47)

Clusters

#Cluster into 3 and get cluster assignments
set.seed(123) #1 2 3 clusters for me
kmeans_clusters <- kmeans(trajectory.curvatures.wide, centers = 3)
clusters <- kmeans_clusters$cluster

kmeans.cluster.results <- data.frame(label = labels, cluster = as.factor(clusters))

 ggplot() + 
   geom_brain(data = kmeans.cluster.results, atlas = glasser, mapping = aes(fill = cluster), colour=I("gray25"), size=I(.08)) + theme_classic() +
 scale_fill_manual(values = c(depth_colorbar[2], depth_colorbar[6], depth_colorbar[3]), na.value = "white") +  
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75"))

## merging atlas and data by 'label'

## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry cluster
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <fct>  
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY 2

## Warning: Some data not merged. Check for spelling mistakes in:
##          label cluster
## 396 lh_L_10pp       2

## merging atlas and data by 'label'

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Clusters_corticalmap.pdf", device = "pdf", dpi = 500, width = 4, height = 2)

## merging atlas and data by 'label'

## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry cluster
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <fct>  
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY 2

## Warning: Some data not merged. Check for spelling mistakes in:
##          label cluster
## 396 lh_L_10pp       2

## merging atlas and data by 'label'

Neurosynth Decoding

neurosynth.terms.clusters <- merge(neurosynth.terms, kmeans.cluster.results, by = c("label"))

clusters.neurosynth <- ldply(neurosynth.termlist, function(t){
  decode.df <- neurosynth.terms.clusters %>% group_by(cluster) %>% do(term.mean = mean(.[[t]])) %>% unnest(term.mean)
  decode.df$term <- t
  return(decode.df)})

clusters.neurosynth %>% filter(cluster == 1) %>% slice_max(n = 5, order_by = term.mean) %>%
ggplot(., aes(x = term.mean, y = term)) +
  geom_segment(aes(y = reorder(term, term.mean), yend = term, x = 0, xend =  
                     term.mean), color = "gray75", linewidth = .5) +
  geom_point(size = 2, color = depth_colorbar[2]) +
  xlab("\nNeurosynth Term Score (z)") +
  ylab("") +
  theme_classic() +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 8, family = "Arial", color = c("black")),
        axis.text.y = element_text(angle = 30, hjust = 1),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

clusters.neurosynth %>% filter(cluster == 3) %>% slice_max(n = 5, order_by = term.mean) %>%
ggplot(., aes(x = term.mean, y = term)) +
  geom_segment(aes(y = reorder(term, term.mean), yend = term, x = 0, xend =  
                     term.mean), color = "gray75", linewidth = .5) +
  geom_point(size = 2, color = depth_colorbar[3]) +
  xlab("\nNeurosynth Term Score (z)") +
  ylab("") +
  theme_classic() +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 8, family = "Arial", color = c("black")),
        axis.text.y = element_text(angle = 35, hjust = 1),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

clusters.neurosynth %>% filter(cluster == 2) %>% slice_max(n = 5, order_by = term.mean) %>%
ggplot(., aes(x = term.mean, y = term)) +
  geom_segment(aes(y = reorder(term, term.mean), yend = term, x = 0, xend =  
                     term.mean), color = "gray75", linewidth = .5) +
  geom_point(size = 2, color = depth_colorbar[6]) +
  xlab("\nNeurosynth Term Score (z)") +
  ylab("") +
  scale_x_continuous(breaks = c(0, 0.5, 1)) +
  theme_classic() +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 8, family = "Arial", color = c("black")),
        axis.text.y = element_text(angle = 30, hjust = 1),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

S-A Axis ranks

kmeans.cluster.results <- merge(kmeans.cluster.results, S.A.axis, by = "label", sort = F)
kmeans.cluster.results %>% group_by(cluster) %>% do(mean.SA = mean(.$SA.axis)) %>% unnest(cols = "mean.SA")

## # A tibble: 3 × 2
##   cluster mean.SA
##   <fct>     <dbl>
## 1 1          148.
## 2 2          278.
## 3 3          241.

Average cluster curvature across depths

cluster.curvatures <- merge(trajectory.curvatures, kmeans.cluster.results, by = c("label")) #curvature for clusters

#Compute average curvature across all regions in a cluster for each cortical depth
cluster.curvatures <- cluster.curvatures %>% group_by(cluster, depth) %>% do(average.curvature = mean(.$curvature)) %>% unnest(cols = average.curvature)
cluster.curvatures$average.curvature.z <- scale(cluster.curvatures$average.curvature, scale = T, center = F)[,1] #scale curvature for ease of plotting. scale divides each data point by the SD
cluster.curvatures$depth <- factor(cluster.curvatures$depth, ordered = T, levels = ordered_depths)

cluster.curvatures %>% 
   ggplot(aes(x = average.curvature.z, y = depth, group = cluster, color = cluster)) +
  geom_point(size = 1.25) +
  geom_path(size = .5) +
  theme_classic() +
  scale_color_manual(values = c(depth_colorbar[2], depth_colorbar[6], depth_colorbar[3]), na.value = "white") +
  xlab("\nCurvature (scaled)") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Depthplot_curvature.pdf", device = "pdf", dpi = 500, width = 2.6, height = 1.6)

Average cluster age of maturation across depths

#calculate average age of maturation at each depth, in each cluster
regional.maturation.alldepths <- merge(regional.maturation.alldepths, kmeans.cluster.results, by = "label", sort = F)
cluster.maturation <-  regional.maturation.alldepths %>% group_by(cluster, depth) %>% do(average.maturation = mean(.$smooth.increase.offset, na.rm = T)) %>% unnest(cols = c("average.maturation"))
cluster.maturation$depth <- factor(cluster.maturation$depth, ordered = T, levels = ordered_depths)

cluster.maturation$cluster.order[cluster.maturation$cluster == 1] <- "a"

## Warning: Unknown or uninitialised column: `cluster.order`.

cluster.maturation$cluster.order[cluster.maturation$cluster == 3] <- "b"
cluster.maturation$cluster.order[cluster.maturation$cluster == 2] <- "c"

cluster.maturation %>% 
  ggplot(aes(y = average.maturation, ymin = 20, x = depth, ymax = average.maturation, group = cluster.order, color = cluster.order)) +
  geom_linerange(position = position_dodge(1), linewidth = 1) +
  scale_y_continuous(breaks = c(20, 22, 24, 26, 28, 30, 32)) +
  ylab("\nAge of Maturation") +
  xlab("Cortical Depth\n") +
  coord_flip() +
  theme_classic() +
  scale_color_manual(values = c(depth_colorbar[2], depth_colorbar[3], depth_colorbar[6]), na.value = "white") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Depthplot_maturation.pdf", device = "pdf", dpi = 500, width = 2.95, height = 1.6)

Curvature variation

trajectory.curvatures$curvature.scaled <- scale(trajectory.curvatures$curvature, scale = T, center = F)
SG.curvature <- trajectory.curvatures %>% filter(depth == "depth_1" | depth == "depth_2" | depth == "depth_3") %>% group_by(label) %>% do(SG.curvature = mean(.$curvature.scaled)) %>% unnest(cols = c("SG.curvature")) #average superficial depths curvature
IG.curvature <- trajectory.curvatures %>% filter(depth == "depth_4" | depth == "depth_5" | depth == "depth_6" | depth == "depth_7") %>% group_by(label) %>% do(IG.curvature = mean(.$curvature.scaled)) %>% unnest(cols = c("IG.curvature")) #average deep depths curaature

SGIG.curvature <- merge(SG.curvature, IG.curvature, by = "label")
SGIG.curvature <- SGIG.curvature %>% mutate(curvature.difference = IG.curvature - SG.curvature) #difference between the two

plot.colorbar <- c("#edf6ff", "#c2d4ed", "#809dc4", "#345c91", "#0e4669")

ggplot() + 
  geom_brain(data = SGIG.curvature, atlas = glasser, mapping = aes(fill = curvature.difference), colour=I("gray25"), size=I(.08)) + 
  theme_classic() + 
  scale_fill_gradientn(colors = plot.colorbar, limits = c(0.4, 1), oob = squish, na.value = "white") +
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75"))

## merging atlas and data by 'label'

## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label          geometry SG.curvature
##   <chr> <chr> <chr> <chr>  <chr> <chr>    <MULTIPOLYGON>        <dbl>
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10…          EMPTY     0.000345
## # ℹ 2 more variables: IG.curvature <dbl>, curvature.difference <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label SG.curvature IG.curvature curvature.difference
## 396 lh_L_10pp 0.0003452231    0.8103573            0.8100121

## merging atlas and data by 'label'

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure3/Curvaturevariation_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.1, height = 2)

## merging atlas and data by 'label'

## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label          geometry SG.curvature
##   <chr> <chr> <chr> <chr>  <chr> <chr>    <MULTIPOLYGON>        <dbl>
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10…          EMPTY     0.000345
## # ℹ 2 more variables: IG.curvature <dbl>, curvature.difference <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label SG.curvature IG.curvature curvature.difference
## 396 lh_L_10pp 0.0003452231    0.8103573            0.8100121

## merging atlas and data by 'label'

Whole-Brain Association Cortex Exhibits Protracted Increases in Superficial Cortex Myelin

Delayed maturation at superficial and all depths

primary.early.sensoryregions <- c("L_V1_ROI", "R_V1_ROI", "L_V2_ROI", "R_V2_ROI", "L_V3_ROI", "R_V3_ROI", "L_V3A_ROI", "R_V3A_ROI", "L_V3B_ROI", "R_V3B_ROI", "L_4_ROI", "R_4_ROI", "L_MT_ROI", "R_MT_ROI", "L_4_ROI", "R_4_ROI", "L_3a_ROI", "R_3a_ROI", "L_3b_ROI", "R_3b_ROI", "L_1_ROI", "R_1_ROI", "L_2_ROI", "R_2_ROI", "L_6d_ROI", "R_6d_ROI", "L_6a_ROI", "R_6a_ROI", "L_FEF_ROI", "R_FEF_ROI", "L_6v_ROI", "R_6v_ROI", "L_6r_ROI", "R_6r_ROI", "L_PEF_ROI", "R_PEF_ROI", "L_A1_ROI", "R_A1_ROI", "L_LBelt_ROI", "R_LBelt_ROI", "L_MBelt_ROI", "R_MBelt_ROI", "L_PBelt_ROI", "R_PBelt_ROI", "L_RI_ROI", "R_RI_ROI") #primary + early sensory regions and primary motor + premotor regions, for outlining in brain plot

glasser.earlysensory <- glasser.regions
glasser.earlysensory$sensory <- ifelse(glasser.earlysensory$orig_parcelname %in% primary.early.sensoryregions, "yes", "no")

#Identify late maturing regions (first, significant increase past age 32)
latematuring <- lapply(gam.statistics.allregions, function(depth){
  depth <- depth %>% filter(smooth.increase.offset > 32) 
  return(depth)
})

#Regions that show a significant increase in R1 past age 32 in both depths 1 and 2
superficial.latematuring <- Reduce(intersect, lapply(latematuring[1:2], `[[`, "orig_parcelname")) %>% as.data.frame %>% set_names("orig_parcelname") %>% merge(glasser.regions) %>% mutate(late.maturation = "yes")
superficial.latematuring <- superficial.latematuring[!(superficial.latematuring$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]

#Plot em
ggplot() +
   geom_brain(data = superficial.latematuring, atlas = glasser,  mapping = aes(fill = late.maturation), colour=I("gray50"), size=I(.06), hemi = "left") + theme_classic() + scale_fill_manual(values = depth_colorbar[3], na.value = "white") + 
   theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
   new_scale_fill() +
   geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex), colour=I(alpha("gray50", 0)), size=I(.06),  hemi = "left") +
   scale_fill_manual(values = c(alpha("white", 0), "gray75")) +
    new_scale_fill() +
   geom_brain(data = glasser.earlysensory, atlas = glasser, mapping = aes(fill = sensory, colour = sensory),  size=I(.12), hemi = "left") +
   scale_color_manual(values = c(alpha("white", 0), "black")) +
   scale_fill_manual(values = c(alpha("white", 0), c(alpha("white", 0))), na.value = "gray75")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure5/Superficial_latematuring_brainplot.pdf", device = "pdf", dpi = 500, width = 2.5, height = 2.5)

#Identify late maturing regions (second, significant increase past age 28)
latematuring <- lapply(gam.statistics.allregions, function(depth){
  depth <- depth %>% filter(smooth.increase.offset > 27) 
  return(depth)
})

#Regions that show a significant increase in R1 past age 32 in both depths 1 and 2
alldepths.latematuring <- Reduce(intersect, lapply(latematuring[1:7], `[[`, "orig_parcelname")) %>% as.data.frame %>% set_names("orig_parcelname") %>% merge(glasser.regions) %>% mutate(late.maturation = "yes")
alldepths.latematuring <- alldepths.latematuring[!(alldepths.latematuring$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]

#Plot em
ggplot() +
   geom_brain(data = alldepths.latematuring, atlas = glasser,  mapping = aes(fill = late.maturation), colour=I("gray50"), size=I(.06), hemi = "right") + theme_classic() + scale_fill_manual(values = depth_colorbar[5], na.value = "white") + 
   theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
   new_scale_fill() +
   geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex), colour=I(alpha("gray50", 0)), size=I(.06),  hemi = "right") +
   scale_fill_manual(values = c(alpha("white", 0), "gray75"))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure5/Alldepths_latematuring_brainplot.pdf", device = "pdf", dpi = 500,width = 2.4, height = 2.4)

Depth-dependent trajectories in exemplar non-frontal regions

gam.smoothestimates.allregions.long <- do.call(rbind, gam.smoothestimates.allregions)
gam.smoothestimates.allregions.long <- gam.smoothestimates.allregions.long %>% mutate(depth = as.factor(substr(row.names(gam.smoothestimates.allregions.long), 1, 7))) #add depth as a column
gam.smoothestimates.allregions.long$depth <- factor(gam.smoothestimates.allregions.long$depth, ordered = T, levels = ordered_depths)

#Function to plot age smooth functions in exemplar regions
plot.depth.smooths <- function(region, y_ticks, line_color){

  smooth.plot <- gam.smoothestimates.allregions.long %>% filter(orig_parcelname == region) %>%
  ggplot(., aes(x = age, y = est, group = depth)) +
  geom_line(linewidth = .3, color = line_color) +
  theme_classic() +
  xlab("\nAge (years)") +
  ylab("R1 Trajectory (zero-centered)\n") +
  scale_y_continuous(breaks = y_ticks) +
  scale_x_continuous(breaks = c(10, 20, 30)) +
  scale_color_manual(values = depth_colorbar) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(smooth.plot)
}

plot.depth.smooths(region = "L_V2_ROI", y_ticks =  c(-0.02, -0.01, 0, 0.01), line_color = "black")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure5/V2_depthsmooths.pdf", device = "pdf", dpi = 500, width = 1.75, height = 1.75)

plot.depth.smooths(region = "L_TE1a_ROI", y_ticks =  c(-0.02, -0.01, 0, 0.01), line_color = depth_colorbar[3])

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure5/TE1a_depthsmooths.pdf", device = "pdf", dpi = 500, width = 1.75, height = 1.75)

plot.depth.smooths(region = "R_7m_ROI", y_ticks = c(-0.015, 0, 0.015), line_color = depth_colorbar[5])

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure5/7m_depthsmooths.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.75)

Charting Depth-Dependent Maturation of Cortical Myeloarchitecture

Valerie Jill Sydnor

Read in Data

Myelin Matures Heterochronously and Hierarchically Across Prefrontal Cortical Depths

Frontal lobe developmental trajectory at each depth

Depth-speficic plots

Multi-depth plots

Number of significant regional effects at each depth

Regional derivatives at each depth

Depth-wise derivative plot

Depth-wise derivative brain plots

Depth correlation matrix

Coefficient of variation

S-A axis correlation

Regional age of maturation at each depth

Depth-wise age of maturation plot

Depth-wise age of maturation brain plots

Depth correlation matrix

Coefficient of variation

S-A axis correlation

Depth-dependent Differences in Myelin Maturational Profiles Vary Across Prefrontal Cortical Zones

Depth-dependent trajectories in exemplar frontal regions

Developmental trajectory curvatures

Clusters

Curvature variation

Whole-Brain Association Cortex Exhibits Protracted Increases in Superficial Cortex Myelin

Delayed maturation at superficial and all depths

Depth-dependent trajectories in exemplar non-frontal regions